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Abstract 

We calculate the strange quark content of the nucleon (iV|ss|iV) in 2 + 1 -flavor lattice QCD. 
Chirally symmetric overlap fermion formulation is used to avoid the contamination from up and 
down quark contents due to an operator mixing between strange and light scalar operators, ss and 
uu + dd. At a lattice spacing a = 0.112(1) fm, we perform calculations at four values of degenerate 
up and down quark masses m u d, which cover a range of the pion mass M n ~ 300-540 MeV. We 
employ two different methods to calculate (N\ss\N). One is a direct method where we calculate 
(iV|ss|./V) by directly inserting the ss operator. The other is an indirect method where (iV|ss|iV) 
is extracted from a derivative of the nucleon mass in terms of the strange quark mass. With 
these two methods we obtain consistent results for (iV|ss|JV) with each other. Our best estimate 
fr s = m s (N\ss\N)/MN = 0.009(15) s t a t(16) S ys is in good agreement with our previous studies in 
two- flavor QCD. 
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I. INTRODUCTION 



The bulk of the nucleon mass Mn is produced by dynamically broken chiral symmetry 
in the vacuum of Quantum Chromodynamics (QCD). This should happen even in the limit 
of vanishing up and down (current) quark masses. Yet, there are also contributions from 
non-zero bare masses of up, down and strange quarks, that are given by a matrix element 
m q (N\qq\N) of a scalar operator qq made of quark field q with mass m q evaluated on the 
nucleon state \N). This quantity is of fundamental importance to characterize the nucleon 
structure. More recently, this quantity, especially that of strange quark, is attracting further 
interest as it determines the cross section of possible dark matter particles to hit the nucleus 
and thus to determine the sensitivity of dark matter search experiments (see, for instance, 

Q). 

The fraction of nucleon mass made of non-vanishing quark masses is conveniently 
parametrized as 

_ m g (N\qq\N) 
lTq M N U 

The light quark contents fx {u d} can be related to the nN sigma term a n N, which is deter- 
mined from experimental data of the ttN scattering amplitude. Evaluation of the strange 
quark content fx 3 is more involved. One uses a n N and a phenomenological estimate of the 
flavor SU(3) violation parameter <jq = m u a(N\uu + dd — 2ss\N), where m u d is (degener- 
ate) up and down quark mass. Recent experimental data a^N = 64(7) MeV {2 1 and <jq = 
36(7) MeV obtained from heavy baryon chiral perturbation theory (HBChPT) {3] led to fx s 
= 0.41(9). This large value appeared to be puzzling, as it suggests that the strange quark 
plays major role to construct nucleon. Early lattice calculations {4-6] also suggested such 
large value. 



In our previous studies 



BE 



8j, we carried out non-perturbative calculations of (N\ss\N) 
in two-flavor QCD, where up and down quarks are assumed to be degenerate. In Ref. [7J, 
(iV|ss|iV) is indirectly estimated from the m s dependence of through the Feynman- 
Hellmann theorem 

Mh m*\n). (2) 



dm. 



is 



on the other 



We refer to this method as the spectrum method in this paper. In Ref. 
hand, (JV|ss|iV) is extracted directly from a disconnected three-point function of the nucleon 
(Fig. [Q. Since we use a ratio of the three- and two-point functions (see fl27j) in Sec. IIIip 



FIG. 1: Nucleoli three-point function used to determine (iV|ss|iV). Solid lines represent quark 
propagators that are dressed by gluons and sea quarks. Connected three lines form the nucleon 
propagator, whereas the disconnected quark loop arising from the strange scalar operator ss. 



to improve the accuracy of (iV|ss|iV), this method is referred to as the ratio method in the 
following. These two studies consistently yielded fx e < 0.05 which is significantly smaller 
than the phenomenological estimate. We note that recent lattice studies (9 J/7 1 also favor 
small strange quark content. 

In this paper, we extend our previous studies to 2 + 1 -flavor QCD. This is a necessary 
step towards a realistic calculation of (iV|ss|iV), since effects of dynamical strange quarks are 
difficult to estimate analytically. In addition, we can eliminate a subtlety in the spectrum 
method when used for two-flavor QCD. Namely, since this theory does not have strange sea 
quark, we estimated (iV|ss|iV) as a derivative in terms of up and down sea quark mass at sea 
{jn u d,sea) and valence (m u d,vai) quark masses set to the physical strange quark mass m SiP h ys 

dM N 



(N\ss\N) 



(3) 

^ lid, sea ^i[d,val ,phys 



9fTLud,sea. 

assuming that (iV|ss|iV) mildly depends on the quark masses. This assumption is eliminated 
in the present work. 
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An important advantage of this work over the previous lattice calculations 
is that chiral symmetry is preserved by employing the overlap quark action 



9|4l9| 



2l|. Con- 



ventional Wilson-type fermions, which explicitly violate chiral symmetry at finite lattice 
spacings, induce a mixing of scalar operators between ss and uu + dd [8|. The nucleon 
three-point function in Fig. [T] then receives a contribution from a connected diagram with 
the uu + dd operator through the renormalization of ss. The connected contribution is larger 
than the disconnected one typically by an order of magnitude, and a subtraction of such a 
large contamination gives rise to a substantial uncertainty in (iV|ss|iV) This serious 

problem is entirely avoided in our work using the chiral lattice fermion formulation. 



This paper is organized as follows. We describe our simulation setup to generate gauge 
ensembles and to calculate relevant nucleon correlators in Sec. [Til The strange quark content 
is extracted through the ratio and spectrum methods at simulated quark masses in Sec. IIHI 
We then extrapolate these results to the physical point in Sec. IIVI Our conclusions are given 
in Sec. |V] Our preliminary reports of this work are found in Refs. 
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II. SIMULATION METHOD 



A. Gauge ensembles 



We simulate QCD with degenerate up and down quarks and heavier strange quarks. 



Chiral symmetry is exactly preserved by employing the overlap quark action 
Dirac operator is given by 
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2lj. Its 



/ 1 1 1 \ / TTl \ 

D(m) = [m + —J + [m - — J 75 sgn[H w (m )), 



(4) 



where m is the quark mass and Hw = ^bDw is the Hermitian Wilson-Dirac operator. The 



mass parameter o: 
has good locality 
proposed in Ref. 



H\y is chosen as mo = — 1.6 so that the overlap-Dirac operator D(m) 
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. For gauge fields, we use the Iwasaki action [25[ with a modification 
261 ] . This leads to an extra Boltzmann factor detfifj^]/ detfi?^ + fi 2 ] 
(/i = 0.2) which does not change the continuum limit of the theory but remarkably reduce 
the computational cost to calculate sgn[if^] in @ by suppressing (near-)zero modes of 
Hw- This Boltzmann factor prohibits tunnelings among different topological sectors, and 
we simulate only trivial topological sector in this study. The effect of fixing topology is 



suppressed by inverse power of the lattice volume 
1 % level, in our previous studies 
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27] and turned out to be small, typically 



291 ] . This small effect can be safely neglected with our 



statistical accuracy for baryon observables. 

Our gauge ensembles are generated at a gauge coupling (3 = 2.3, where the lattice spacing 
is determined as a = 0.112(1) fm using the Q baryon mass as input. On a N^xN t = 16 3 x 48 
lattice, we simulate two values of the degenerate up and down quark masses m u d = 0.035 
and 0.050, and two strange quark masses m s = 0.080 and 0.100. Their physical values 
m ud, P hys = 0.0029 and m SjP hys — 0.081 are fixed by using M w and M K as inputs [29]. Note also 
that we quote bare values in lattice units for these quark masses. We push our simulations to 
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0.015 


0.025 


0.035 


0.050 


m s 


0.080 0.080 0.100 


0.080 0.080 0.100 


0.080 0.100 


0.080 0.100 


L/a 


24 16 16 


24 16 16 


16 16 


16 16 



TABLE I: Summary of parameters used in the lattice simulation. 

two smaller m U dS, 0.015 and 0.025, on a larger lattice 24 3 x 48 at a single value of m s = 0.080, 
which is very close to m sphys . 

Four values of m u d cover a range of the pion mass ~ 300-540 MeV. The spatial 
extent L is chosen to satisfy a condition M V L > 4 to control finite volume effects. We carry 
out additional simulations at two smallest m u dS on the smaller lattice size 16 3 x48 to directly 
examine the finite volume effects. Our simulation parameters are summarized in Table HI 

The statistical samples at each simulation point (m u d,m s , L) consist of 2,500 hybrid 
Monte Carlo trajectories, out of which we use 500 and 50 to calculate the correlation func- 
tions in the spectrum and ratio methods, respectively. We employ the jackknife method 
with a bin size of 50 trajectories to estimate statistical errors of the nucleon correlators and 
any quantities determined from them. 

On these gauge ensembles, we calculate the two-point nucleon correlation function using 
an interpolating operator N = e ahc {u^C^dh)u c with C = 7472- After taking contractions, 
we obtain 

<C 2pt (y,t,At)> = --L £ J2^ a ' b ' c ' 

s r=(i± 74 )/2 x 

(|tr4r(D- 1 (m))->r4r(D- 1 (m)) 66 '(C'75)(( J D- 1 (m)) cc ') r (C75)] 

+tr s [r( J D- 1 (m))-'(C75)(( J D- 1 (m)) cc ') T (C75)p- 1 (m)) bb '] V (5) 

where the trace "tr s " is over spinor indices and (■ • • ) represents a Monte Carlo average. Here, 
the quark propagators D~ l (m) propagate from (y, t) to (x, t + At). In order to improve 
statistical accuracy, C 2pt is averaged over two choices of the projector T = (1 ±7 4 )/2, which 
correspond to the forward and backward propagating nucleons, respectively. Here and in 
the following, for F — (1 — j±)/2, At is taken as —At. 

We also calculate the three-point function with a scalar operator on the lattice defined 
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as 

llat ' D 



2mo / 

to respect chiral symmetry in the continuum limit. 



B. All-to-all propagator 

As shown in Fig. [TJ the three-point function Cs pt on a given gauge configuration can be 
decomposed into two pieces. Namely, we can write Cs pt as 

(C 3pt {y,t,At,At s )} = (C 2pt (y,t,At)S la \t + At s )), (7) 

where C 2pt (y,t, At) is the two-point function and 



S^(t + At s ) = ^^^Tt(D-\m))(z,z)\ M0= ^M. ~ (Tr(D-\m))(z, z)\ Zo=t+Ats ^ } 



(8) 



is the scalar quark loop calculated on this configuration. The trace "Tr" is over both spinor 
and color indices. The nucleon piece C 2p t can be calculated by using the conventional 
"point-to-all" quark propagator /} _1 (x,x'), the source point of which (x') has to be fixed 
to a certain lattice site. The calculation of the quark-loop pieces S' lat is computationally 
more demanding, as it involves quark loops starting from arbitrary lattice sites (z,t + At s ). 



We therefore employ the "all-to-all" quark propagator 30|, |3l| that contains the quark 
propagating from any lattice site to any site. 

Let us consider a decomposition of the quark propagator to the contribution from low- 
lying eigenmodes of the Dirac operator D(m) and that from the remaining modes 

D-\m) = {D-\m)} Xow + {D' 1 (m)}^. (9) 

It is expected that the low-mode contribution {D~ 1 (m)}i ow dominates for low-energy ob- 
servables in QCD including (N\ss\N). We calculate it exactly as 

{D _1 (m)}i ow (s, y) = y~] v k (x)v k {y)\ (10) 

ti A»(m) 

where Afc(m) and Vk(x) are the k-th lowest eigenvalue and its associated eigenvector of D(m), 
and iV e is the number of the low-lying modes prepared for this calculation. 

The small contribution from the remaining high-modes is calculated stochastically by the 



noise method 



32]. We generate a single complex Z 2 noise vector r)(x) for each configuration 



and split it into Nd — 3 x 4 x N t /2 vectors rj^ d \x) (d = l,...,Nd), which have nonzero 
elements only for a single combination of color and spinor indices on two consecutive time 
slices. For each "split" noise vector r]^ d \ we solve a linear equation 

{D(m)^}(x) = (V hish V (d) ){x) (d=l,...,N d ), (11) 
where Phi g h = 1 — 7\>w and V\ ow is the projector to the subspace spanned by the low-modes 

Vi ow (x, y) = ^2 v k( x ) v k(yY- ( 12 ) 

k=l 

The high-mode contribution is then estimated as 

N d 

{D-\m)} high (x,y) = ^^(^W- (13) 

We calculate the low- and high-mode contributions to S lat as 

S lat (t + At.) = Sl:l(t + At s ) + Slf gh (t + At,), (14) 

with 

s{:i {high) (t+At s ) = ^X^^^HWigh)^^)!^^, (is) 

s z 

where the subtraction of the VEV is assumed here though it is not written explicitly for 
notational simplicity. 

C. Low- mode averaging (LMA) 

The low-lying modes of D(m) are also useful to precisely calculate the nucleon piece C^pt 
in both C 2p t and C 3pt . By applying (jUJ), we can decompose C 2v t into the following eight 
contributions 

/~i r^lll I /~i\lh I rilhl < rihll , rilhh i r^hlh i s~ihhl , rihhh (~\(\\ 

^2pt — <^2pt ' °2pt ' °2pt ' °2pt ' °2pt "+" °2pt "+* 2pt "+* Ly 2pt • 1 10 J 

Here, C^. is constructed only from {D _1 (m)} low . For C^pt; {-^ _1 ( m )}iow is used for two of 
the valence quark propagators and {-D _1 (?7i)}high for the remaining one. The other combina- 
tions are understood in a similar manner. In principle, we can use the all-to-all propagator, 



S 



(fit)]) and ({TBI) , to calculate these contributions. These quantities however decay exponen- 
tially with a large nucleon mass as the temporal separation At increases. At large 
separations, the high- mode contributions, such as C^, are not sufficiently precise with 
{Z} _1 (m)}hi g h evaluated using only single noise sample for each configuration. 



We therefore use the low-mode averaging (LMA) technique [33, |34j in this study. The 
low- mode part of the all-to-all propagator (TTOT) is used to calculate C^L, which dominantly 
contributes to the nucleon correlators C*2 P t and C^pt- We then take average of C^iy, t, At) 
over the location of the nucleon source operator (y, t) to largely reduce its statistical fluc- 
tuation. 

The remaining and small contributions {C^, C^/ 1 } are calculated using the point-to- 
all quark propagator after projecting by V\ ow and l — V\ ow for I and h pieces, respectively. 
We improve the statistical signal of these contributions by averaging over (y,t). In order 
to reduce the computational cost of the re- calculation of the point-to-all propagators, these 
contributions are averaged over a limited set of (y,t) compared to that for C^L- 

The sets of the source point as well as the number of the low-modes N e are chosen 
differently for our calculations with the ratio and spectrum methods, because the latter 
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29|. We 



uses C 2p t calculated in the course of our study of the light meson spectrum 
summarize our choices for these two methods in the following subsections. 

D. Setup for the ratio method 

We use N e = 160 and 240 low-lying modes on the 16 3 x48 and 24 3 x48 lattices, respectively, 
to calculate low-mode contribution S^lit + At) and C l 2p t (y,t, At) in the ratio method. As 
mentioned above, the latter is averaged over 16 spatial points 

y G {(0, 0, 0), (0, 0, N s /2), (0, N s /2, N s /2), (N s /2, N s /2, N s /2), (N s /4, N s /4, N s /4), 
(AV 4, A,/ 4, 3 Ay 4), (N 8 /4, 3N s /4, 3A s /4), (3A s /4, 3A s /4, 3JV./4), 
and their permutations} (17) 

at each time slice t. Averaging over more points does not help to further reduce the statistical 
fluctuation of C 2p t and Cs p t because of the correlation among C2p t (y, t, At) at different spatial 
points y's. We average {C^pt; C^//} over four time slices t = 0, 12, 24 and 36 with the 
spatial location y kept fixed. 
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At heaviest m u d (=0.050), we slightly modify the setup of LMA to calculate Cs pt . With 
()14p and f fl6|) . C^ vt on a given configuration can be rewritten as 

olat I J /~illh I I f~ihhh \ rdat 

3pt — <^2pt D low "T \ ( ~ y 2pt "T " " ' ~r <-^2pt J D low 

"T°2pt D high "T t°2pt "'*""' °2pt J °high- 1 10 J 

In general, C3 P t is largely dominated by the first term. The other three terms have a totally 
small contribution due to a significant cancellation among their statistical fluctuations. As 
m u d increases, the first term less dominates C^ pt . We observe that the cancellation among 
the remaining terms becomes more significant, which is however spoiled by LMA, namely 
averaging only the third term over more choices of (y,t), leading to a larger error of C% v t- 
We therefore apply LMA only for the first term at the largest m uc i- 

The above setup of LMA typically leads to a factor of 4 (7) reduction of the statistical 
error of C*2 P t (C3 P t) at our simulated values of m U( i and m s . 

In our previous study in two-flavor QCD [8], we observe that smearing both nucleon 
source and sink operators is crucial to identify the ground state contribution to Cz pt . We 
employ the Gaussian smearing 

«CM) = E { ( 1 + W H Y) *(y,f), H ^y = IX y -; + W> ( 19 ) 

for both C*2pt and Cz pt after fixing the gauge. The parameters u = 20 and iV = 400 are 
chosen by inspecting the plateau of the effective mass of C^pt- 



E. Setup for the spectrum method 

For the spectrum method, we use C 2pt calculated in the course of our study of the 
light meson spectrum 



28|, [29|. The low-mode contribution C^ t (y, t, At) is calculated using 
iV e = 160 (80) low-modes on the 16 3 x48 (24 3 x48) lattice, and is averaged over the time-slice 
t with the spatial source point y kept fixed. We use an exponential smearing 

CS(x, *) = E exp[-0.4|r|]g(x + r, t) (20) 

r 

only for the nucleon source operator. The spatial extent of this smeared operator is roughly 
equal to that of (fT9|) used for the ratio method. We observe that the onset of the plateau in 
the effective mass is consistent with that of ffl9l) within the statistical error. 
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In order to evaluate the derivative DMn I dm K in ([2]), we study the m s dependence of M, 



by utilizing the reweighting technique 



35 



N 



36J. Our Monte Carlo data at the strange quark 



mass m s are used to estimate the two-point function at a slightly shifted strange quark mass 



m s as 



(C2pt)m' s 



(C 2pt w(m' 8l ffl 8 )), 



(21) 



where (• • ■) ms represents the Monte Carlo average at m s , and w is the reweighting factor 
for a given configuration 



w m Q , m, 



w(m' s , m s ) 



det 



Dim. 



(22) 



(w(m' s ,m s )) ms ' 

Similarly to 5* lat and C*2 P t, w can be decomposed into contributions from low- and high- 
modes 



w(m' s ,m s ) = wi ow (m' s ,m s )w high (m' s ,m s ), 

Dim'; 



Wlow(high)(ra s ,m s J 



det 



low(high) 



Dim, 



low(high) 



(23) 
(24) 



We exactly calculate w\ ow using the low-lying eigenvalues, whereas Whigh is estimated by a 
stochastic estimator for its square 

N r 



w high {m s ,m s ) 



Nr r=l 



(n-l)Phigh^r 



(25) 



Here Q = D(m s )\D(m' s )~ i y D^m'^^Dirris), and {Ci; CiV r } i s a se ^ °f pseudo-fermion 
fields whose elements are generated with the Gaussian probability. 

An important practical issue is how many pseudo-fermion fields are needed to reliably 
estimate ifhigh- Since u>hi g h is a product of 12N^N t — N e eigenvalues, it largely deviates 
from unity unless m s ~m' s . We observe, however, that it has small statistical fluctuation, 
after taking the ratio w(m' s , m s ) = w{m' s , m s ) / '(w(m' s , m s )) ms . Consequently, the normalized 
reweighting factor w is essentially controlled by the low-mode contribution u>i ow . We there- 
fore do not need large number of the pseudo-fermion fields to estimate Whigh as demonstrated 
in Fig. H 

In this study, we reweight C 2p t at m s = 0.080 to 20 different values 

m ' s = 0.0600, 0.0650, 0.0700, 0.0725, 0.0750, 0.0775, 0.0780, 0.0785, 0.0790, 0.0795, 

0.0805, 0.0810, 0.0815, 0.0820, 0.0825, 0.0850, 0.0875, 0.0900, 0.0950, 0.1000.(26) 
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°1bm 1500 2000 

trajectory 

FIG. 2: Monte Carlo history of reweighting factor w(m' s ,m s ) to shift the strange quark mass 
from m' s = 0.080 to m s = 0.075 at m u d = 0.050. Different lines show data calculated with different 
numbers of the pseudo-fermion fields N r . 

We shift these values by +0.020 when we reweight C*2 P t at m s = 0.100. These values roughly 
cover a region m' s G [m s — 25 MeV, m s + 25 MeV] , where the low- mode dominance of w is 
confirmed. We set N r = 5 in the whole region of m' s . 



III. RESULTS AT THE SIMULATED QUARK MASSES 

In the following subsections, we present our results for (NlO^lN) obtained at simulated 
quark masses by using the ratio and spectrum methods. Note that (N\Og\N) represents 
the bare value on the lattice, and results for the renormalization group invariant parameter 
Jt s will be given in the next section. 



A. Ratio method 

We extract (N\Of\N) from the ratio of C 3pt (At, At s ) and C 2pt (At) 

R(At, At.) = C3 * (A *;f*' } —T > (N\Of\N), (27) 

where At is the temporal interval between the nucleon source and sink. The scalar quark loop 
S lat is set on the time-slice apart from the nucleon source by At s . Note that Cs p t(At, At s ) 
and C2 P t(A£) are calculated using LMA and, hence, we suppress the coordinates of the 
nucleon source, namely (y,t) in (jSJ) and (jjTJ). 



12 




FIG. 3: Approximated ratio R(At, At s ) at m s = 0.080 as a function of At s . We plot results 
obtained at different values of m u d in the four panels. The vertical dashed lines show the locations 
of the nucleon source and sink operators. 



The ratio R may receive contamination from excited states of the nucleon when the 
temporal separation At is not sufficiently large and/or the scalar operator is too close to the 
nucleon operators (At s ~ or At). We therefore need to identify a plateau of R(At, At s ), 
where the excited state contamination is negligible. To this end, we consider the same ratio 
but approximated by taking only the low-mode contribution for the quark loop S la,t in 
( !T4|) . This approximated ratio, which we denote by R in the following, is useful to identify 
the plateau of R, because i) R is well dominated by the low-mode approximation R, and ii) R 
is free from a large noise due to the stochastic method to estimate S^ gh , which obscures the 
excited state contamination. We refer the reader to Ref. [8] for a more detailed discussion. 

Figure [3] shows R(At, At s ) with a fixed value of At = 13 as a function of At s . We obtain 
nonzero signal for R(At, At s ), which do not show significant At s dependence at At s ~ At/2. 
It implies that the scalar operator is sufficiently far from nucleon operators. 
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FIG. 4: Results of the constant fit to R(At, At s ) as a function of At. Four panels show results 
obtained at different values of m u d and at m s = 0.080. 



We then carry out a constant fit to the ratio without the approximation R(At, At s ) using 
a fit range of At s = [5, At — 5} for each At. As plotted in Fig. HJ the fit results are stable at 
At > 12 where our data are well dominated by the ground state contribution. 

^From these observations on the At s and At dependences, we determine (N\Og\N) by a 
simultaneous constant fit to R(At, At s ) with fit ranges of At s = [5, At— 5] and At = [12, 23]. 
The numerical results are listed in Table [Til We also test a fitting form taking account of 
the first excited state with a slightly wider fit range of At s . This fit yields (N 1 O l g"\ N) in 
good agreement with those from the constant fit, because the excited state contribution is 
small as expected from the mild At dependence of R(At, At s ). 

We repeat the same analysis at two smallest m^'s but on the smaller volume 16 3 x48. 
The numerical results are listed in Table [TTTI The difference in (N\Og\N) between the two 
volumes are well below our statistical accuracy suggesting that finite volume effects (FVEs) 
can be neglected within the statistical error. We therefore use the numerical results in 



14 



m ud 


0.015 


0.025 


0.035 0.035 


0.050 0.050 


m s 


0.080 


0.080 


0.080 0.100 


0.080 0.100 


(N\O l ^\N) 


0.09(7) 


0.15(6) 


0.28(12) 0.14(10) 


0.15(14) 0.20(18) 



TABLE II: Strange quark content (N\Og\N) calculated in the ratio method. 





0.015 0.015 


0.025 0.025 


m s 


0.080 0.100 


0.080 0.100 


(N\Of\N) 


0.34(24) 0.29(32) 


0.21(16) 0.02(9) 



TABLE III: Same as Table M but for m ud = 0.015 and 0.025 on the smaller volume 16 3 x48. 
Table HT1 in the chiral extrapolation to determine (NlO^lN) at physical quark masses. 



B. Spectrum method 

In the spectrum method, we evaluate (NlO^lN) from the m s dependence of the nucleon 
mass Mn- Figure [5] shows examples of the nucleon effective mass obtained at m s = 0.080. 
By a single exponential fit C*2 P t(At) oc e _MjvAi , we determine Mjv with an accuracy of 2% 
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FIG. 5: Nucleon effective masses at m s = 0.080. Left and right panels show results on the 16 3 x 48 
and 24 3 x 48 lattices, respectively. Horizontal lines show Mjy obtained from a single exponential 
fit to C 2pt (At). 
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FIG. 6: Nucleon masses Mjv as a function of m s . Left and right panels show Mn obtained by 
reweighting from that at m s = 0.080 and 0.100, respectively. We plot results on the 16 3 x 48 and 
24 3 x 48 lattices, by open and filled symbols. Filled squares are slightly shifted in the horizontal 
direction for clarity. 
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FIG. 7: Fitted results for {N\0^\N} as a function of the width of the fitting range Am s . Left and 
right panels show results at (m u d, m s ) = (0.015, 0.080) and (m u d, m s ) = (0.050, 0.080), respectively. 

(0.8%) at our smallest (largest) m ud on 24 3 x 48 (16 3 x 48). 

The FVE in M/y at m ud = 0.025 is not statistically significant: it is only la (3%) effect. 
We expect similarly small effect at heavier m Mrf 's. The magnitude of the FVE at m ud = 0.015 
is difficult to estimate due to a large statistical error of M N on the smaller volume 16 3 x 48. 
We note that the FVE at m ud = 0.015 on the larger volume 24 3 x 48 is estimated as 0.7% 
from SU(2) heavy baryon chiral perturbation theory (HBChPT) at one-loop. In addition, it 
is plausible that the FVE has a mild dependence on m s leading to small effect in (iVjO^liV). 
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0.015 


0.025 


0.035 0.035 


0.050 0.050 


m s 


0.080 


0.080 


0.080 0.100 


0.080 0.100 


(N\O l ^\N) 


-0.16(35) 


0.35(13) 


0.31(15) 0.16(12) 


0.42(10) 0.22(10) 



TABLE IV: Strange quark content (N\Og\N) obtained from the spectrum method. 





0.015 


0.015 


0.025 0.025 


m s 


0.080 


0.100 


0.080 0.100 


(N\of\N) 


0.68(52) 


-0.60(55) 


0.34(11) 0.12(18) 



TABLE V: Same as Table HV] but at m ud = 0.015 and 0.025 on the smaller volume 16 3 x48. 

As explained in the previous section, we calculate at shifted values of m s by ex- 
ploiting the reweighting technique. Results are plotted as a function of m s in Fig. El We 
successfully reweight our data to m s ± 0.02 (±25 MeV). Namely, the reweighting does not 
largely increase the statistical error of M^. This is because i) the reweighting factor w, 
is accurately calculated with the small number of the noise samples, as discussed in the 
previous section, and ii) resulting values are typically 0(1) as plotted in Fig. [2j 

We extract the slope dM^ / dm s by fitting in the region of [m s — Am s ,m s + Am s ] 
with Am s = 0.010 to a linear form 

M N = d+ (N\Of\N)m s . (28) 

The numerical results are summarized in Table. IIV1 Figure [7] shows that the fitted result for 
( A^O^I N) is stable against the choice of the fitting range Am s as expected from the mild 
m s dependence of shown in Fig. |6j We also confirm that adding higher order terms to 
( 128]) does not change (NlO^lN) significantly. 

In order to directly check FVEs to (JVlOg^iV), we repeat the analysis at two lightest m U d 
but on the smaller volume. A comparison with results listed in Table |V] suggests that FVE 
is not significant with our statistical accuracy, which is consistent with our observation in 
the ratio method. 
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FIG. 8: Strange quark content (JV|C?g t |iV") as a function of m u d (given in the lattice unit). Left 
and right panels show the results from the spectrum and ratio methods, respectively. 

C. Comparison between two methods 

Figure [8] compares (NlO^lN) obtained from the spectrum and ratio methods. We ob- 
serve a good agreement between the two methods. The same figure also shows that FVEs 
in (NlO^lN) are not significant at the two smallest m,,/s as already mentioned in the 
previous subsections. These observations suggest that systematics of our determinations at 
given quark masses (m u d, rn s ) is not substantial. 

With our simulation setup, the accuracy at two heaviest m u dS are comparable between 
the two methods, while the ratio method provides a more accurate determination at lighter 
m u dS. This is mainly because i) we use the better setup of LMA for the ratio method and 
ii) the volume size is increased at these m u dS. For instance, we average Cj^y, t, At) at the 
16 choices of the spatial location y, while y is kept fixed in the spectrum method. Our data 
at the 16 choices listed in (|T71) have less correlation among them on a larger volume and 
hence LMA works better. 

IV. CHIRAL EXTRAPOLATION 

In Fig. [9X we plot (N\Og\N) obtained from the two methods as a function of m U( i. Note 
that our data cover a region of M n ~ 300 - 540 MeV, and our lighter strange quark mass 
m s = 0.080 is already close to the physical mass m s phys = 0.081. The figure shows that 
(NlO^lN) has a very mild dependence on both m u d and m s , which has also been observed 
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FIG. 9: Linear extrapolations of (N\Og\N) obtained from ratio (left panel) and spectrum (right 
panel) methods. Solid and dashed lines show the fit lines at m s = m StP ^y S and 0.100, respectively. 
We omit the fit line at m s = 0.080, which is not indistinguishable from that at m s = WJ S) phys m the 
scale of the figure. Star symbols represent (NlO^lN) extrapolated to the physical quark masses. 
Square symbols are slightly shifted in the horizontal direction for clarity. 



X 2 /d.o.f. d.o.f. 


C Ci 


ci, s (N\Of\N) 


spectrum method 0.54 3 


0.90(47) 5.7(5.1) 


-9.6(5.4) 0.15(0.19) 


ratio method 0.38 3 


0.22(43) 3.9(4.1) 


-2.2(5.7) 0.058(0.101) 



TABLE VI: Numerical results of linear chiral extrapolation, 
in our previous study in two- flavor QCD [8|. Our data are well described by a linear fit 

llat 



(N\Of\N) = c + c 1M m ud + c 1)S m s 



(29) 



as plotted in the same figure. Numerical results of the fit are summarized in Table I VI I We 
also confirm that (N\Og\N) at the physical quark masses does not change significantly by 
excluding the data at the largest m u d from the fit and/or by including higher order terms 
in (EH]). 

We also test a fitting form based on SU(3) HBChPT to parametrize the observed quark 
mass dependence of (iV|ss|iV). One- loop chiral expansion of 37| and the Feynman- 
Hellmann theorem ([2]) gives an expression of (iV| |iV) 

(N\Of\N) = -c s -B \^C NNK M K + 2C NNn M„| + c 2<K M 2 K + c^M 2 ^ (30) 

where contributions of the decuplet baryons are ignored. In this analysis, we approximate 
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FIG. 10: Chiral fits of {N\O l f\N) using 5(7(3) HBChPT (30]). Left and right panels show fits to 



(NlO^lN) obtained from the ratio and spectrum methods, respectively. 



X 2 /d.o.f. d.o.f. — c s 


c 2 ,K [GeV- 2 ] 


c 2>v [GeV- 2 ] 


(N\O l ^\N) 


spectrum method 1.2 3 8.7(5) 


23(4) 


-6.1(3.6) 


0.33(19) 


ratio method 0.75 3 7.9(4) 


21(3) 


-2.8(3.6) 


0.16(10) 



TABLE VII: Numerical results of chiral fit using SU(3) HBChPT (1301) . 



the higher order corrections by the 0(M? K analytic terms. Within this approximation, 
we can use the leading order expressions M\ = B{m ud + m s ) and M 2 = 2B{m ud + 2m s )/3 
for the meson masses. The coefficients C^nk and C^n v of the 0(M{K,rj\) terms are written 
as 

1 (5.D 2 — QDF + 9F 2 ) 1 (D — 3F) 2 



C 



NNK 



C 



8vr/ 2 3 ' ~ NN " 8ti P 6 

The axial couplings are fixed to a phenomeno logical estimate D = 0.81 and F = 0.47 [38j 
in this analysis. The low-energy constants in mesonic ChPT, / and B, are set to our lattice 
estimate determined from the meson spectrum and decay constants [29]. 

The fit using (130!) is shown in Fig.[10l Since (iV|ss|iV) depends mildly on m u d through the 
strange meson masses Msk }7] \ up to one-loop order of HBChPT, the mild m u d dependence 
of our data can be fitted to (130!) reasonably well. However, numerical results summarized in 
Table IVHl suggest a large difference of (NlO^lN) in the SU (3) chiral limit between the linear 
and HBChPT fits (cf. —c s in Table [VTT1 and cq in Table IVT1) . This is because ( 1301) predicts a 
large O(Mk) contribution to (iVlOg^iV) at m SiP h ys with the phenomenological estimate of 
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FIG. 11: Contribution to (N\O l ^\N) in the chiral expansion (|30p . Each contribution are calculated 
at the physical strange quark mass. 

D and F. Then, the fit reproduces our small values of (NlO^lN) by a large cancellation 
among chiral corrections at different orders. Consequently, the HBChPT expansion exhibits 
a poor convergence as shown in Fig. [TTJ A similarly poor convergence of HBChPT has been 
observed in our study in two-flavor QCD [8(. These observations suggest that, at least for 
(iV|Og ,t |JV), the SU(3) chiral expansion up to 0(M? K y) could be applicable only to lattice 
data at much smaller values of m s . 

In this study, therefore, we determine (N\Og\N) from the linear fit (!29j) and use the 
HBChPT fit only to estimate the systematic uncertainty of the chiral extrapolation. We 
obtain (N\O l ^\N) = 0.15(19) (18) from the spectrum method and 0.06(10)(10) from the ratio 
method. The first and second errors represent the statistical and systematic ones. In this 
study, the ratio method provides a statistically better determination of (N\ss\N). This 
is partly because we employ a better setup of LMA for the ratio method as mentioned in 
subsections III Dl and IIII Cl A nucleon operator with a better overlap with the nucleon ground 
state also improves the accuracy of the spectrum method. 

As discussed in Sec. IHH we expect that the FVE on our larger volume is small. The 
discretization effect is estimated as 0((aA) 2 ) ~ 9% from a simple order counting using A = 
500 MeV. These systematic errors are well below our statistical accuracy and, hence, ignored 
in the following discussions. We also note that exact chiral symmetry in our simulation, 
forbids the mixing with the light quark contents (N\uu + dd\N) 
to introduce a large uncertainty in (iV|ss|./V) 18]. 
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181 ] , which turned out 



The bare matrix element (N\Og\N) is converted to the renormalization invariant pa- 
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N f = 2 + 1 and 2+1+1 QCD |sHig|. We convert m s (N\ss\N) or (N\ss\N) in 



Q, 14, Q 



to St. 



using the experimental value of Mn and m s obtained in [391 ] . 



rameter 



h 



m s {N\Of\N) 



M 



N 



0.023(29)(28) (spectrum method), 
0.009(15)(16) (ratio method), 



(32) 



where we use experimental value of Mm. In Fig. [T^J we compare our results of St s with 
our previous estimate in N f = 2 QCD Q, fl. All of our studies give consistent results for 



St s ■ As confirmed in Fig. [131 sea strange quark loops have small effects to a renormalization 
invariant quantity m s (N\O l g t \N) leading to the good agreement between Nf = 2 and 2+1 
QCD. As mentioned in the Introduction, our previous study using the spectrum method 
in Nf = 2 QCD |7J estimated (-/V|ss|iV) from the derivative dM N / dm ud ^ seil with m nrfj { seaiVa i} 
sending to m s phys . This turns out to be a reasonable estimate of (iV|ss|iV) because of the 
very mild dependence of (JV|ss|iV) on ms^s} shown in Fig. [H] as well as the small effect of 
dynamical strange quark loops in Fig. [TBI 

Figure [T2l also compares our results with recent studies in Nf = 2 + 1 and 2+1+1 lattice 
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FIG. 13: A comparison of a re-normalization invariant quantity m s (N\ss\N) between Nf = 2 
(Refs. [8], open symbols) and 2+1 (this study, filled symbols) QCD. We plot data from the ratio 
method as a function of Ml. 



QCD |9Hl6j|. All these studies favor small strange quark content fr s < 0.1. Strictly speaking, 
the results of Ref. [l0| appears to be slightly higher (2.5 a) than our best estimate, that is 
fx s in Nf = 2 + 1 QCD from the ratio method. Recently, same authors present improved 



estimates in Nf = 2 + 1 and 2 + 1 + 1 QCD 15(. These results also indicate a slightly large 
value of fx a ~ 0.06. Given the large statistical errors, however, the difference is not very 
significant. 

Compared to these lattice estimates, the phenomenological studies {2, 3] predict a rather 
large estimate 0.41(9) based on HBChPT up to quadratic order in the quark masses. The 
poor convergence of our chiral fit based on the same effective theory suggests that its con- 
vergence at physical quark masses should be carefully examined. 



V. CONCLUSION 



In this paper, we calculate the strange quark content of the nucleon in 2 + 1 -flavor lattice 
QCD. Two determinations using the ratio and spectrum methods as well as our previous 
studies in two-flavor QCD consistently favor a small strange quark content fx B < 0.05. In 
contrast, phenomenological studies based on HBChPT have led to a rather larger value 
0.1-0.7 which is, however, unexpectedly large as a content of sea quarks of a single flavor. 
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In this study, we utilize several simulation techniques to precisely determine the small 
effect due to disconnected quark loops. In the spectrum method, we can successfully shift 
m s by ±25 MeV by using the reweighting technique. It would be interesting to study isospin 
breaking effects, such as the proton and neutron mass difference, by using this technique, 
and its applicability on larger lattice volumes should be studied. 

The ratio method requires precise calculation of the nucleon disconnected three-point 
function, which is technically very challenging. The low-lying modes of the Dirac operator 
turned out to be very helpful: we employ the LMA technique to calculate the nucleon 
propagator and the all-to-all quark propagator for the disconnected quark loops. These 
techniques, in principle, can be applied to other baryon observables. For instance, it is 
interesting to extend this study to the strange quark spin content of the nucleon. Precise 
knowledge of this quantity is important to constrain the parameter space of SUSY models 
through spin-dependent scattering cross section of the neutralino- nucleon scattering 40 1. 
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